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Abstract 



Ewald summation and physically equivalent methods such as particle-mesh 
Ewald, kubic-harmonic expansions, or Lekner sums are commonly used to 
calculate long-range electrostatic interactions in computer simulations of po- 
lar and charged substances. The calculation of pressures in such systems is 
investigated. We find that the virial and thermodynamic pressures differ be- 
cause of the explicit volume dependence of the effective, resummed Ewald 
potential. The thermodynamic pressure, obtained from the volume derivative 
of the Helmholtz free energy, can be expressed easily for both ionic and rigid 
molecular systems. For a system of rigid molecules, the electrostatic energy 
and the forces at the atom positions are required, both of which are read- 
ily available in molecular dynamics codes. We then calculate the virial and 
thermodynamic pressures for the extended simple point charge (SPC/E) wa- 
ter model at standard conditions. We find that the thermodynamic pressure 
exhibits considerably less system size dependence than the virial pressure. 
From an analysis of the cross correlation between the virial and thermody- 
namic pressure, we conclude that the thermodynamic pressure should be used 
to drive volume fluctuations in constant-pressure simulations. 
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I. INTRODUCTION 



Pressure is one of the fundamental thermodynamic variables. The calculation of pressures 
in fluid systems using computer simulations is generally considered to be a routine task. 
However, difficulties arise in the presence of long-range interactions. Here, we investigate 
the calculation of pressures in computer simulations of charged and polar systems, where 
the long-range Coulomb interactions are commonly treated with Ewald lattice summationB 
or physically equivalent methods like particle-mesh Ewald,i kubic-harmonic expansions,^ 
or Lekner sums.S^ A mechanistic definition of the pressure leads to the standard virial 
expression. A thermodynamic definition of the pressure is based on the volume dependence 
of the Helmholtz free energy. When the Coulomb interactions are resummed by using, e.g., 
the Ewald method, the resulting effective pair interactions depend explicitly on the volume. 
In addition, self interactions are present that also depend on the volume. As a consequence, 
the virial and thermodynamic pressures are not identical for finite Coulomb systems, even 
though the two pressures are expected to converge in the thermodynamic limit. 

The paper is organized as follows: in section ^ we derive expressions for the virial and 
thermodynamic pressures. In section |T|, we study the pressure in systems with long-range 
Coulomb interactions. For the thermodynamic pressure, we derive a simple formula that can 
be readily implemented in standard molecular dynamics or Monte Carlo codes. In sections 

and 0, we study the system size dependence of the virial and thermodynamic pressures 
for the extended simple point charge (SPC/E) water mo delS under standard conditions. 



II. VIRIAL AND THERMODYNAMIC PRESSURES 
A. Virial pressure 

The pressure p can be calculated from a mechanistic prescription equating the exterior 
and interior forces on the container. This leads to the virial expression for the pressure in 
an atomic system,! 
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Pv = P^bT + — 




(1) 



where p is the number density of particles; fee is Boltzmann's constant; T is the temperature 
{kBT = l3~^); and V is the volume. The sum extends over the scalar product betwen particle 
positions rj and forces Fj exerted on particle i due to other particles in the system. (. . .) 
denotes a canonical average. For computer simulations under periodic boundary conditions 
with pair forces, it is convenient to rewrite the virial equation in a form that makes the 
translational invariance explicit: 



where r^j = r^- — and Fjj = —dv{Yij) / dvi = dv{Yij) / dvij is the pair force exerted on 
particle i by particle j, derived from a pair potential f (r); and the sum is over all pairs of 
particles in the system. 

For a system of rigid polyatomic molecules i,j with interaction sites a and /?, one obtains 
an analogous formula when the forces Yiajp between molecular sites are projected onto a 
vector between the "centers" of the two molecules (e.g., the center of mass). 



Here, Fjj is the net force between two rigid molecules, summed over molecular sites a and 
j3. Otherwise, the constraint forces maintaining the rigidity of the molecules have to be 
included explicitly in Eq. @. 



The thermodynamic expression for the pressure is derived from the relation between the 
pressure pr, the Helmholtz free energy F, and the volume V , 




(2) 




(3) 



B. Thermodynamic pressure 
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(4) 



T,N 



The statistical-mechanical relation between the free energy F and the partition function 



Qn{V,T) for identical classical particles in a canonical ensemble is: 



QN{V,T) = e 



-I3F _ 



1 



// 



(5) 



where h is Planck's constant; H = K + U is the Hamiltonian; and the integration extends 
over the positions and momenta of all particles. In taking the volume derivative 
Eq. the kinetic energy K is independent of the volume. Transforming the positional 
coordinates into dimensionless form, V~^r'^ , and pulling out a factor from the integral, 
leads to the ideal gas term pk^T for the pressure. The non-ideal contributions are contained 
in the volume dependence of the potential energy U, 



Typically, U does not depend explicitly on the volume. The volume dependence of U then 
arises from the volume scaling of the particle positions. In the absence of an explicit volume 
dependence, we can express dU / dV as 



with dvijdV = Ti/SV. By using Fj = —dU/dVi and combining Eqs. (H) and (0), we find the 
corresponding thermodynamic pressure to be equivalent to the mechanistic pressure Eq. (p. 




(6) 




I 1 



(7) 



III. PRESSURE IN SYSTEMS WITH LONG-RANGE COULOMB 



INTERACTIONS 



A. Thermodynamic pressure in ionic systems 



The identity between the virial and thermodynamic pressures, Eqs. (|I]) and (^, does 
not hold if the potential depends explicitly on the system volume. Such an explicit volume 
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dependence arises in computer simulations of charged and polar substances, when the long- 
range Coulomb interactions are resummed, e.g., by using the Ewald method.§ 

We first split the total potential energy U into a short-range part [/(^'') and a long-range 
Coulomb part U^^\ 

U = U^""'^ + f/(") . (8) 

The pressure p^^^^ corresponding to U^^"^^ contains the ideal-gas term and the contributions 
from short-range pair interactions, 

P^^^^ = P^BT--L^^Fg,.r..' (9) 

where the short-range forces 'Pfajp those derived from the short-range part U^""'^ of the 
potential energy. Note that the virial and thermodynamic expressions for p'^^^^ are equivalent, 
and therefore the subscripts "F' or "T" have been omitted in Eq. (^). 

Next, we consider the pressure arising from the potential energy f/(^) of Ion g-range 
Coulomb interactions. In Ewald lattice summation, the charges in a periodically repli- 
cated simulation box interact with an effective potential. That potential is obtained from a 
summation over all periodic images. In addition, a self interaction arises from interactions 
with a particle's own images. This leads to a Coulomb energy U^"^^ for a system of partial 
charges qia at positions rjQ,: 

a, 13 



•ifi 



a<l3 



i a 



I ^iaifS I 
1 



(10) 



The first sum is the intermolecular contribution; the second and third sums are the in- 
tramolecular contributions, with the self interactions contained in the third sum. (f{r) is the 
effective, resummed Coulomb potential, with a Fourier representation:!^ 
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V'(r) = ^ E P"--' , (11) 

where the k sum extends over the reciprocal lattice corresponding to the lattice vectors n of 
periodically replicated simulation boxes. In a cubic lattice of length L = V^^^, we have n = 
L k), and k = 27rL"^ k), where i, j, and k are integers. It is numerically convenient 
to partly transform ip{r) into real space, leading to its Ewald lattice sum representation, 

^ erfc(r,|r+ n|) ^ 4. _ . 

n |r + n| ^ Kk"^ Kr/^ 

?7 is a convergence parameter chosen to accelerate numerical convergence. The value of (f{r) 
is independent of ?7,0 

The self-interactions in [/^^^ are given by the interactions of a unit point charge with its 
periodic images, subtracting the bare self interaction, (p{r) — l/|r|, with the appropriate 
limit taken for r ^ 0. For a given box shape, (f{r) scales with the box volume V as 

^(r) = \/-V3^*(r*) , (14a) 
^ = --^(r) , (14b) 

where star superscripts denote volume-independent quantities. This follows from Eq. ([IlD 
with volume scaling r ~ V^^^ and k ~ V^^^^. The same scaling is true trivially for the direct 
1/1 r I interactions. For an ionic system of point charges without bond constraints, Eqs. 
and (p!^ immediately lead to an expression for the thermodynamic pressure in terms of the 
Coulomb energy 

PT = P^-^ + ^ . (15) 

Equation (p^ ) gives the well-known relation between the pressure and energy of an ionic 
system, for which the Coulomb energy is a homogeneous function of degree —1 in the 
coordinates.0 



B. Thermodynamic pressure in systems of rigid polyatomic molecules 



For a system of rigid molecules, we find the following volume scaling: 



dip{Yiajp) _ 1 



dV 
_d 1_ 

dV Ir,™ 



iai/3 1 



3V 
1 

W 



1 d 1 



HaifS 



= 



(16a) 
(16b) 

(16c) 



where Vij is the distance vector between two molecule centers; dj^ = — rj is the vector 
from the center to site a; and diajp = dj/j — djQ,. Equation (|16c|) follows from the volume 
independence of the intramolecular distance vector rj^j^ = djQ,j^. Combining Eqs. (p!oD and 
(p!6|), we find for the volume derivative of the Coulomb energy: 



5f/(c) f/(c) 



dV 



3V 3V 



EE 

»a a, 13 



d 



iaj/3 



iajfS) 



+ EE 

a<f3 



d 



1 



iaifi ) 



I ^iaifS I 



d-iajfj 



Half} 



(17) 



We can simplify dU^'^^ /dV further by expressing it in terms of the intermolecular forces 



''iajj^ exerted by site j(3 onto site ia, 



(inter) _ 



d 



(18) 



and the intramolecular forces F^-^*^'^'', 



p (intra) 

iaifi 



d 



dr 



iaifS 



-qiaqip 



iai/3 ) 



^iai/3 1 _ 



p (intra) 

il3ia 



(19) 



This leads to 



dV 



1 



Jjic) + ^ ^ F.^J^J-* ■ diajf3 + E E ^iaip"^ ■ dim/3 | , 



(20) 



i,j a, (3 

i<3 



i a, 13 
a<l3 



The sums over pairs of sites ia and j(3 can be rewritten as a single sum over all sites. This is 
possible because the distances diajp = djf3 — dia are intramolecular and are continuous when 
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a particle crosses the box boundary (i.e., diajp does not change when the periodic images of 
the particles i or j are used) . This results in 

5f/(c) 1 



dV 3V 

where F^^**^'^'' is the net intermolecular force on site ia. 



(21) 



FSr^ = EEFSS^ (22) 



and F^"*'^'^'' is the net intramolecular force on site ia, 



-p(mtra} _ ^ -p(mtra) ^23") 

/3 

It is advantageous to add the intra and intermolecular forces because in common Ewald-sum 
implementations the Fourier term already contains the sum of both inter and intramolecular 
contributions which are thus not easily separated. We define the net Coulomb force F^ on 
site ia as the sum of the inter and intramolecular forces, 

Tp(c) _ Tp(inter) -p(intra) _ (2A) 
^ ia ^ ia ' ia p, ' V / 

^^ia 

We then find for the thermodynamic pressure of a system of rigid molecules: 

PT=P^-~^ + ^ ((f/(^^) - (EEfS ■ d..^) . (25) 

Thus the presence of intramolecular constraints in rigid polyatomic molecules resulted in a 
force term to be subtracted from the pressure of the purely ionic system, Eq. (|1^). Note 
that the forces F^^ in Eq. (^) are derived from the Coulomb energy U''^^ alone. Additional 
ideal-gas and short-range contributions to the pressure are reflected in p^^'^\ 



C. Tin- foil boundary conditions and reaction field correction 



The infinite Ewald lattice is implicitly embedded in a conducting medium with dielectric 
constant e^f = cxd, corresponding to "tin-foil" boundary conditions. This is the appropriate 
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choice for a conducting medium. However, for a polar substance it can be advantageous to 
use a reaction-field dielectric constant e^f similar to that of the bulk medium.lll The resulting 
correction to the Coulomb energy U^^\ey-{ = oo) is0 

271 



(26) 



(2erf + l)y' 

M is the instantaneous dipole moment of the simulation volume arising from the dipole 
moments irij of individual molecules, 



(27) 



For rigid molecules the rrij do not change with volume. The reaction-field correction, 
Eq. (p6D, thus scales as V~'^, 

5f/(rf) f/{rf) 



dV V ■ 

The forces derived from the reaction-field correction are 

^(rf) _ ^[/(■•f) _ Anq,^ 



M . 



(28) 



(29) 



dr,^ (2erf + 1)V 

By using Eq. (pT]), we can express the sum of reaction- field forces projected onto the in 
tramolecular distance vectors in terms of the reaction-field energy 



i:FS?.d,„ = -2t/(^^). 

Accordingly, the volume derivative of the reaction-field energy can be written as 



(30) 



5f/{rf) 



(31) 



dV 3V 

The correction Eq. (|26| ) for a finite reaction-field dielectric constant e^f then leads to an 
expression for the thermodynamic pressure analogous to Eq. (p5|). 



3V 



i a 



Pt = P 

Here, the forces F^^,'^" are derived from the Coulomb energy 

f/(c) 



(32) 



p(c,erf) 



5f/(c) I 



erf) 



and contain the reaction field contribution F-^, defined in Eq. (^ 



dri, 

(rf) 



(33) 
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D. General considerations 



We emphasize the simphcity of the pressure expressions Eqs. (p5[) and (|32[) for systems 
of rigid molecules. The Coulomb energy contribution ijj'^'^'^) /ZV is analogous to that of 
the corresponding ionic system, Eq. (0), corrected for the presence of constraint forces. A 
more or less equivalent expression for the pressure in Coulombic systems treated with Ewald 
summation was derived before by Smith,0 and similarly by Boulougouris et al^. However, 
in those derivations the volume derivative was carried out explicitly for the Ewald energy. 
Also, the derivations start from an approximate expression for the Ewald energy that does 
not include the full real-space lattice sum and self terms. Therefore, the derivations did not 
arrive at a closed expression and the simplicity of the results given here was masked. 

Expressions for the pressure tensor P for Ewald summation were derived previously by 
Nose and Klein,Ei and Heyes,0 as discussed by Alejandre et a/., El as well as by others.EHnl 
However, the tensor character does not lend itself easily to a compact notation for the 
bulk pressure p = Tr(P) in a homogeneous system. Equations (|1^), (p5D, and (|32D have 
the advantage of being independent of the specific method used to evaluate the energies 
and forces. All that is needed is the total Coulomb energy and forces at all sites that are 
consistent with that energy. This is what molecular dynamics codes will normally produce at 
no additional cost. The Coulomb interactions can then be evaluated by using conventional 
Ewald sums,0il particle-mesh Ewald,i kubic-harmonic expansion!, or Lekner sums.i^ For 
approximate Coulomb energy calculations such as reaction-fieldii'ii or generalized reaction- 
field methods,!! Eqs. (|^), (^Sf), and (^) suggest an evaluation of the pressure that is 
formally consistent with that of Ewald sums and physically equivalent methods. 



IV. COMPUTER SIMULATIONS 

To investigate the quantitative differences between the virial and thermodynamic pres- 
sures, we study a model of water at standard conditions (298 K temperature, 997.07 kg m~^ 
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mass density corresponding to a number density of p = 33.33 nm ^). We use the SPC/E 
model of water,! formed by a Lennard- Jones center on the oxygens, 

^LJ(-) = ^2--e^ (34) 

where A = 0.3428^2 kJ nm^Vmol and B = 0.371226 ^^j nmVmol. In addition, the SPC/E 
model carries three partial charges. The hydrogen and oxygen sites carry charges qh = 
0.4238e and qo = — 2g//, respectively, where e is the elementary charge. The oxygen- 
hydrogen bond length is 0.1 nm, the hydrogen-oxygen-hydrogen bond angle is cos"^(— 1/3) ~ 
109.47 deg. 

We use Metropolis Monte Carlo simulations for the canonical sampling,^! where the 
translational and rotational move widths are chosen to give an acceptance rate of about 
40 per cent. Ewald summation is used for the electrostatic interactions, with rj = 5.6/ L 
where L = V^^^ is the length of the cubic box. A spherical cutoff of L/2 is used for the 
real space interactions (charge and Lennard- Jones). The real-space potentials are shifted 
by a constant, such that they are zero at the cutoff. The Fourier space sum is truncated at 

< 38(27r/L)^, leading to 2 x 510 k vectors being considered. A reaction-field dielectric 
constant of erf = 65 has been used in all simulations. Standard finite-size corrections were 
applied to the Lennard- Jones contributions to pressure and potential energy.il 

System sizes of = 16, 32, 64, 128, 256, and 512 water molecules are studied. Starting 
from random configurations, these systems have been equilibrated for at least 250 000 Monte 
Carlo passes. (One pass corresponds to one attempted move for each of the A^ particles.) 
In the production runs, the energy as well as the virial and thermodynamic pressures are 
calculated every tenth pass. 

V. RESULTS FOR SPC/E WATER 

Table | contains the simulation characteristics, as well as results for the virial and ther- 
modynamic pressures. The thermodynamic pressure is calculated using Eq. (^). The virial 
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pressure is calculated using Eq. @, where the pair forces are derived from the total potential 
energy U = U^'^'^ + f/("). Also included in Table Pl are results for the average potential en- 
ergy per particle. (To compare with the experimental energy, one has to add a polarization 
correction of 5.22 kJ/mol.@) Errors are obtained from a block analysis,ii plotting calculated 
standard deviations of the mean as a function of the number of blocks used. The reported 
error is then the plateau value reached in the limit of long blocks with typically more than 
about 10 000 Monte Carlo passes. 

The system size dependence of the virial and thermodynamic pressure is shown in Fig. |l|. 
From Table | and Fig. |l], we find that pv and pt converge to the same value of about —5 
MPa (1 MPa = 10 bar) for large system sizes, with a statistical error of about 2 MPa. This 
convergence is expected as the thermodynamic and virial pressure should be identical in the 
thermodynamic limit. However, the thermodynamic pressure exhibits a considerably weaker 
system size dependence than the virial pressure. The thermodynamic pressure for as few as 
64 SPC/E water molecules is in agreement with large system sizes. The virial pressure, on 
the other hand, scales as roughly l/N"^ for small to intermediate system sizes, with its value 
off by about one order of magnitude for = 64. We emphasize that for typical system sizes 
of > 256, the virial and thermodynamic pressures are identical within statistical errors 
for SPC/E water under standard conditions. 

Figure shows the radial distribution functions of water oxygens and hydrogens, which 
were calculated also in the corners of the cubic simulation box with appropriate weights. We 
find that the the radial distribution functions for A^ > 64 water molecules are practically in- 
distinguishable, whereas the A^ = 16 and A^ = 32 simulations are somewhat more structured 
beyond the first peaks. These shght structural differences could explain the deviations of 
the thermodynamic pressure for those small system sizes. We caution that these are results 
for the specific thermodynamic state (room temperature and standard density) studied here, 
and we expect more pronounced finite-size effects, e.g., for low densities. 

In constant pressure simulations,0@ the box volume is rescaled according to the "in- 
stantaneous pressure" obtained from individual configurations by omitting the canonical 
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average (. . .) in the pressure formulas above. It is therefore important that not only the 
average but also the instantaneous pressure driving the volume fluctuations be correct. As 
measures of discrepancy between the virial and thermodynamic pressures pv and pt, we use 
the correlation coefficient r and the average absolute deviation minus the deviation of the 
averages, A, 

^ _ {{PT - {pt)){pv - {pv))) 

{{PT - {pt)YY'\{Pv - {pv)??" ' 
^ = {\Pt-{pt)-Pv + {Pv)\) , (35b) 

where instantaneous pressures pt and pv are used. Results for r and A are listed in Table |. 
The cross-correlation coefficient r indicates strong correlation, with r values between 0.90 
and 0.997 for = 16 to iV = 512. However, the average absolute deviation A between the 
two pressures is significant even for systems of 512 water molecules, scaling approximately 
as A ~ 1/A^. Therefore, in constant pressure simulations, the use of the thermodynamic 
pressure appears advantageous. 

In an earlier study of pressure effects on the stability of hydrophobic aggregates in 
water we determined the thermodjTiamic pressure of SPC wateS as a function of density 
using Eq. P^). For the temperature and density studied here (T = 298 K, p = 33.33 nm~^), 
we found a pressure of about 37 ± 6 MPa for SPC water. From the density dependence of 
the pressure, we determined a compressibility factor phsTxr ~ 0.06 for SPC water, where 
Xt is the isothermal compressibility. That compressibility factor is in excellent agreement 
with the experimental value of 0.062. 



VI. CONCLUSIONS 

We have derived a simple, compact expression for the Coulomb contribution to the ther- 
modynamic pressure in a system treated with Ewald lattice summation. For a system of 
point ions, we recover the well-known relation between the pressure and potential energy. 
We then derive an expression for the pressure in a system of rigid molecules carrying point 
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charges. The pressure in such a system can be calculated from the total energy and the forces 
at each site alone. This makes the implementation of that pressure formula trivial, because 
both energy and forces are produced routinely in molecular dynamics codes. Moreover, these 
formulas are entirely independent of the particular method used to resum the Coulomb inter- 
actions. Ewald summation, particle-mesh Ewald,li kubic-harmonic expansions,! and Lekner 
sumJi^ can be used readily. For approximate reaction-field methods,&ii expressions for 
the pressure are suggested by analogy. 

We have compared the thermodynamic pressure, obtained from the volume dependence 
of the Helmholtz free energy, with the mechanistic virial pressure. We find that for rigid 
SPC/E water at standard conditions, the two pressures are approximately equal (within 
errors of about 2 MPa) for systems larger than = 256 water molecules. For smaller 
systems, the virial pressure exhibits a pronounced system-size dependence, whereas the 
thermodynamic pressure can be calculated accurately by using as few as 64 SPC/E water 
molecules. 
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FIGURES 

FIG. 1. Pressure of SPC/E water as a function of the inverse number of water molecules, 1/N. 
Cross symbols and dashed lines correspond to the virial pressure pv- Plus symbols and solid lines 
show the results for the thermodynamic pressure pr- The inset highlights results for larger system 
sizes, N > 64. Error bars indicate one standard deviation of the mean, estimated from a block 
error analysis. 

FIG. 2. Radial distribution functions of water atoms. Oxygen-oxygen (top panel), oxy- 
gen-hydrogen (middle panel) , and hydrogen- hydrogen (bottom panel) radial distribution functions 
are shown for different numbers of water molecules. Arrows indicate half the box length, r = L/2, 
for different system sizes. The radial distribution functions were calculated for distances beyond 
L/2 using appropriate weights. 
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TABLES 

TABLE L Characteristics and results of the Monte Carlo simulations of SPC/E water. Statisti- 
cal errors are one estimated standard deviation of the mean. Also included are the cross-correlation 
coefficient r and the absolute deviation A, as defined in Eq. (|35|). 





passes [10'^] 


{U/N) [kJ/mol] 


Pv [MPa] 


Pt [MPa] 
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A [MPa] 
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-28.6 ±5.6 
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